library('mvtnorm')
plot.mixture <- function(model, YY, ...) {
plot(YY, ...)
mmu <- model$mmu
SS <- model$SS
m <- length(mmu)
xx <- seq(min(YY[,1]), max(YY[,1]), by=.1)
yy <- seq(min(YY[,2]), max(YY[,2]), by=.1)
for (ell in 1:m) {
contour(xx, yy,
outer(xx, yy, function(xx, yy) dmvnorm(cbind(xx, yy), mmu[[ell]], SS[[ell]])),
add=T, col=1+ell%%7)
}
}
YY <- as.matrix(faithful)
par(mfrow=c(1,2))
# boxplot(YY$eruptions)
# boxplot(YY$waiting)
par(mfrow=c(1,1))
plot(YY)
llk <- function(model, YY) {
mmu <- model$mmu
SS <- model$SS
pp <- model$pp
m <- length(mmu)
N <- nrow(YY)
PPHI <- matrix(0, N, m)
for (ell in 1:m) {
PPHI[,ell] <- dmvnorm(YY, mean=mmu[[ell]], sigma=SS[[ell]])
}
sum(log(PPHI %*% pp))
}
em.gauss <- function(m, YY, llk.diff.eps=0.1, do.plot=FALSE, verbose=FALSE) {
d <- ncol(YY)
N <- nrow(YY)
phi <- dmvnorm
mmu <- list()
SS <- list()
ww <- list()
zz <- list()
llk.prev <- -Inf
llk.hist <- NULL
# Initialization
# TODO: Split the data into m clusters at random; calculate mmu, SS, pp for each cluster
for (ell in 1:m) {
# pick two data points at random
mmu[[ell]] <- sample(as.data.frame(t(YY)), 1)[[1]]
SS[[ell]] <- cov(YY)
}
pp <- rep(1/m, m)
if (do.plot) {
plot.mixture(list(pp=pp, mmu=mmu, SS=SS), YY, main=paste('n.iter =', 0))
}
n.iter <- 1
repeat {
# Expectation
PPHI <- matrix(0, N, m)
for (ell in 1:m) {
PPHI[,ell] <- phi(YY, mean=mmu[[ell]], sigma=SS[[ell]])
}
for (ell in 1:m) {
zz[[ell]] <- pp[[ell]] * PPHI[,ell] / PPHI %*% pp
}
for (ell in 1:m) {
ww[[ell]] <- zz[[ell]] / sum(zz[[ell]])
}
# Maximization
for (ell in 1:m) {
mmu[[ell]] <- as.vector(t(YY) %*% ww[[ell]])
pp[[ell]] <- sum(zz[[ell]]) / N
# FIXME: vectorize
for (i1 in 1:d) {
for (i2 in 1:d) {
SS[[ell]][i1, i2] <- sum(ww[[ell]] * (YY[, i1] - mmu[[ell]][i1]) * (YY[, i2] - mmu[[ell]][i2]))
}
}
}
llk.cur <- llk(list(pp=pp, mmu=mmu, SS=SS), YY)
llk.diff <- abs(llk.cur - llk.prev)
if (verbose) {
cat('# ', n.iter, ' LLK = ', llk.cur, ' (', llk.diff, ')\n', sep='')
}
llk.hist <- c(llk.hist, llk.cur)
if (do.plot) {
plot.mixture(list(pp=pp, mmu=mmu, SS=SS), YY, main=paste('n.iter =', n.iter))
}
if (llk.diff < llk.diff.eps) {
if (do.plot) {
plot(llk.hist, type='o', col='blue', main='Log likelihood', xlab='n.iter')
}
return(list(pp=pp, mmu=mmu, SS=SS, zz=zz, llk=llk.cur, llk.hist=llk.hist))
} else {
llk.prev <- llk.cur
n.iter <- n.iter + 1
}
}
}
set.seed(2)
model.2 <- em.gauss(2, YY, do.plot=TRUE, verbose=TRUE)
## # 1 LLK = -1266.242 (Inf)
## # 2 LLK = -1235.165 (31.0769)
## # 3 LLK = -1207.979 (27.1865)
## # 4 LLK = -1200.716 (7.263274)
## # 5 LLK = -1196.433 (4.282454)
## # 6 LLK = -1192.152 (4.281588)
## # 7 LLK = -1186.858 (5.294085)
## # 8 LLK = -1179.499 (7.358276)
## # 9 LLK = -1168.903 (10.59576)
## # 10 LLK = -1154.545 (14.3581)
## # 11 LLK = -1140.924 (13.62177)
## # 12 LLK = -1131.893 (9.030789)
## # 13 LLK = -1130.32 (1.572512)
## # 14 LLK = -1130.267 (0.05378514)
model.2$pp
## [1] 0.6438597 0.3561403
model.2$mmu
## [[1]]
## [1] 4.290234 79.974944
##
## [[2]]
## [1] 2.037046 54.485312
model.2$SS
## [[1]]
## eruptions waiting
## eruptions 0.1692466 0.9315412
## waiting 0.9315412 35.9457512
##
## [[2]]
## eruptions waiting
## eruptions 0.06969635 0.440863
## waiting 0.44086299 33.739059
set.seed(1)
model.3 <- em.gauss(3, YY, do.plot=TRUE, verbose=TRUE)
## # 1 LLK = -1286.069 (Inf)
## # 2 LLK = -1285.762 (0.3072853)
## # 3 LLK = -1285.387 (0.3742914)
## # 4 LLK = -1284.823 (0.5646781)
## # 5 LLK = -1283.835 (0.9871156)
## # 6 LLK = -1281.761 (2.074725)
## # 7 LLK = -1276.226 (5.534951)
## # 8 LLK = -1258.127 (18.09884)
## # 9 LLK = -1219.377 (38.75022)
## # 10 LLK = -1199.715 (19.66186)
## # 11 LLK = -1188.31 (11.40481)
## # 12 LLK = -1175.743 (12.56663)
## # 13 LLK = -1161.195 (14.54839)
## # 14 LLK = -1145.295 (15.8997)
## # 15 LLK = -1132.672 (12.62365)
## # 16 LLK = -1127.521 (5.15066)
## # 17 LLK = -1123.44 (4.080583)
## # 18 LLK = -1121.309 (2.131676)
## # 19 LLK = -1120.531 (0.7774306)
## # 20 LLK = -1120.231 (0.3000882)
## # 21 LLK = -1120.087 (0.1446344)
## # 22 LLK = -1119.999 (0.08806301)
Choose the best model out of several runs with the highest LLK.
fit.em.gauss <- function(m, YY, n.tries=10, verbose=FALSE) {
prev.llk <- -Inf
model.best <- NULL
for (n.iter in 1:n.tries) {
if (verbose) {
cat('= Fitting model', n.iter, '...\n')
}
model <- em.gauss(m, YY, verbose=verbose)
if (model$llk > prev.llk) {
model.best <- model
}
}
return(model)
}
set.seed(1)
model <- fit.em.gauss(3, YY, verbose=TRUE)
## = Fitting model 1 ...
## # 1 LLK = -1286.069 (Inf)
## # 2 LLK = -1285.762 (0.3072853)
## # 3 LLK = -1285.387 (0.3742914)
## # 4 LLK = -1284.823 (0.5646781)
## # 5 LLK = -1283.835 (0.9871156)
## # 6 LLK = -1281.761 (2.074725)
## # 7 LLK = -1276.226 (5.534951)
## # 8 LLK = -1258.127 (18.09884)
## # 9 LLK = -1219.377 (38.75022)
## # 10 LLK = -1199.715 (19.66186)
## # 11 LLK = -1188.31 (11.40481)
## # 12 LLK = -1175.743 (12.56663)
## # 13 LLK = -1161.195 (14.54839)
## # 14 LLK = -1145.295 (15.8997)
## # 15 LLK = -1132.672 (12.62365)
## # 16 LLK = -1127.521 (5.15066)
## # 17 LLK = -1123.44 (4.080583)
## # 18 LLK = -1121.309 (2.131676)
## # 19 LLK = -1120.531 (0.7774306)
## # 20 LLK = -1120.231 (0.3000882)
## # 21 LLK = -1120.087 (0.1446344)
## # 22 LLK = -1119.999 (0.08806301)
## = Fitting model 2 ...
## # 1 LLK = -1234.707 (Inf)
## # 2 LLK = -1154.135 (80.57213)
## # 3 LLK = -1128.975 (25.15987)
## # 4 LLK = -1127.784 (1.190921)
## # 5 LLK = -1126.533 (1.251448)
## # 6 LLK = -1125.433 (1.099579)
## # 7 LLK = -1124.573 (0.8599739)
## # 8 LLK = -1123.894 (0.678957)
## # 9 LLK = -1123.329 (0.565005)
## # 10 LLK = -1122.858 (0.4711058)
## # 11 LLK = -1122.491 (0.3669589)
## # 12 LLK = -1122.223 (0.2675441)
## # 13 LLK = -1122.026 (0.1975011)
## # 14 LLK = -1121.868 (0.1579442)
## # 15 LLK = -1121.73 (0.137739)
## # 16 LLK = -1121.603 (0.1276629)
## # 17 LLK = -1121.48 (0.1225137)
## # 18 LLK = -1121.36 (0.1196227)
## # 19 LLK = -1121.243 (0.1176387)
## # 20 LLK = -1121.127 (0.1158781)
## # 21 LLK = -1121.013 (0.1140001)
## # 22 LLK = -1120.901 (0.111843)
## # 23 LLK = -1120.792 (0.1093402)
## # 24 LLK = -1120.685 (0.1064765)
## # 25 LLK = -1120.582 (0.1032649)
## # 26 LLK = -1120.482 (0.09973455)
## = Fitting model 3 ...
## # 1 LLK = -1277.145 (Inf)
## # 2 LLK = -1259.889 (17.25532)
## # 3 LLK = -1218.301 (41.58858)
## # 4 LLK = -1179.961 (38.33938)
## # 5 LLK = -1145.735 (34.2264)
## # 6 LLK = -1124.045 (21.68996)
## # 7 LLK = -1121.814 (2.230887)
## # 8 LLK = -1120.946 (0.8686486)
## # 9 LLK = -1120.557 (0.3884166)
## # 10 LLK = -1120.364 (0.1927393)
## # 11 LLK = -1120.24 (0.1243156)
## # 12 LLK = -1120.142 (0.09763961)
## = Fitting model 4 ...
## # 1 LLK = -1272.919 (Inf)
## # 2 LLK = -1242.463 (30.45639)
## # 3 LLK = -1202.45 (40.01349)
## # 4 LLK = -1183.286 (19.16338)
## # 5 LLK = -1161.098 (22.1881)
## # 6 LLK = -1133.807 (27.29108)
## # 7 LLK = -1125.533 (8.27415)
## # 8 LLK = -1122.386 (3.146323)
## # 9 LLK = -1120.903 (1.483316)
## # 10 LLK = -1120.117 (0.7864931)
## # 11 LLK = -1119.695 (0.42118)
## # 12 LLK = -1119.475 (0.2209002)
## # 13 LLK = -1119.36 (0.1149018)
## # 14 LLK = -1119.299 (0.06034246)
## = Fitting model 5 ...
## # 1 LLK = -1288.577 (Inf)
## # 2 LLK = -1288.068 (0.5084829)
## # 3 LLK = -1287.597 (0.4707323)
## # 4 LLK = -1286.995 (0.6021216)
## # 5 LLK = -1286.026 (0.9691306)
## # 6 LLK = -1284.114 (1.911752)
## # 7 LLK = -1279.376 (4.738229)
## # 8 LLK = -1264.716 (14.66069)
## # 9 LLK = -1228.145 (36.57005)
## # 10 LLK = -1206.73 (21.41557)
## # 11 LLK = -1199.818 (6.911756)
## # 12 LLK = -1194.515 (5.303474)
## # 13 LLK = -1187.251 (7.263575)
## # 14 LLK = -1174.669 (12.58193)
## # 15 LLK = -1153.767 (20.90242)
## # 16 LLK = -1131.888 (21.87843)
## # 17 LLK = -1123.439 (8.449541)
## # 18 LLK = -1120.886 (2.552951)
## # 19 LLK = -1120.271 (0.6144042)
## # 20 LLK = -1120.081 (0.1908042)
## # 21 LLK = -1119.989 (0.09145503)
## = Fitting model 6 ...
## # 1 LLK = -1288.16 (Inf)
## # 2 LLK = -1287.846 (0.3134081)
## # 3 LLK = -1287.449 (0.397249)
## # 4 LLK = -1286.875 (0.5743079)
## # 5 LLK = -1285.943 (0.9320087)
## # 6 LLK = -1284.168 (1.774968)
## # 7 LLK = -1279.969 (4.198858)
## # 8 LLK = -1268.352 (11.61639)
## # 9 LLK = -1248.586 (19.76627)
## # 10 LLK = -1234.311 (14.27505)
## # 11 LLK = -1220.973 (13.33768)
## # 12 LLK = -1209.253 (11.72008)
## # 13 LLK = -1199.989 (9.264551)
## # 14 LLK = -1191.219 (8.769459)
## # 15 LLK = -1178.88 (12.33982)
## # 16 LLK = -1157.449 (21.43039)
## # 17 LLK = -1130.774 (26.67494)
## # 18 LLK = -1121.587 (9.186763)
## # 19 LLK = -1120.036 (1.551722)
## # 20 LLK = -1119.734 (0.3017132)
## # 21 LLK = -1119.621 (0.1131207)
## # 22 LLK = -1119.55 (0.0708383)
## = Fitting model 7 ...
## # 1 LLK = -1287.116 (Inf)
## # 2 LLK = -1286.855 (0.2609048)
## # 3 LLK = -1286.738 (0.1176217)
## # 4 LLK = -1286.673 (0.06533044)
## = Fitting model 8 ...
## # 1 LLK = -1277.105 (Inf)
## # 2 LLK = -1264.75 (12.3548)
## # 3 LLK = -1231.585 (33.16466)
## # 4 LLK = -1180.806 (50.77975)
## # 5 LLK = -1159.923 (20.88298)
## # 6 LLK = -1146.501 (13.42163)
## # 7 LLK = -1132.846 (13.65529)
## # 8 LLK = -1128.071 (4.774265)
## # 9 LLK = -1127.709 (0.3620281)
## # 10 LLK = -1127.572 (0.1376114)
## # 11 LLK = -1127.436 (0.1360656)
## # 12 LLK = -1127.292 (0.143773)
## # 13 LLK = -1127.14 (0.1523581)
## # 14 LLK = -1126.979 (0.1608453)
## # 15 LLK = -1126.81 (0.1688024)
## # 16 LLK = -1126.634 (0.1759478)
## # 17 LLK = -1126.452 (0.1821249)
## # 18 LLK = -1126.265 (0.1872763)
## # 19 LLK = -1126.073 (0.1913882)
## # 20 LLK = -1125.879 (0.1944187)
## # 21 LLK = -1125.683 (0.1962423)
## # 22 LLK = -1125.486 (0.1966388)
## # 23 LLK = -1125.291 (0.1953392)
## # 24 LLK = -1125.098 (0.1921203)
## # 25 LLK = -1124.912 (0.1869126)
## # 26 LLK = -1124.732 (0.179869)
## # 27 LLK = -1124.56 (0.1713491)
## # 28 LLK = -1124.398 (0.1618245)
## # 29 LLK = -1124.247 (0.1517568)
## # 30 LLK = -1124.105 (0.1415104)
## # 31 LLK = -1123.974 (0.1313287)
## # 32 LLK = -1123.853 (0.1213551)
## # 33 LLK = -1123.741 (0.111671)
## # 34 LLK = -1123.639 (0.1023272)
## # 35 LLK = -1123.545 (0.09336253)
## = Fitting model 9 ...
## # 1 LLK = -1233.099 (Inf)
## # 2 LLK = -1186.704 (46.39474)
## # 3 LLK = -1164.808 (21.89594)
## # 4 LLK = -1145.595 (19.21353)
## # 5 LLK = -1134.166 (11.42824)
## # 6 LLK = -1126.752 (7.414515)
## # 7 LLK = -1123.039 (3.712398)
## # 8 LLK = -1121.495 (1.544885)
## # 9 LLK = -1120.837 (0.6574851)
## # 10 LLK = -1120.508 (0.3292761)
## # 11 LLK = -1120.311 (0.1967074)
## # 12 LLK = -1120.178 (0.1327207)
## # 13 LLK = -1120.082 (0.09668214)
## = Fitting model 10 ...
## # 1 LLK = -1246.34 (Inf)
## # 2 LLK = -1201.032 (45.3079)
## # 3 LLK = -1179.25 (21.78218)
## # 4 LLK = -1161.035 (18.21422)
## # 5 LLK = -1139.702 (21.33387)
## # 6 LLK = -1126.924 (12.7777)
## # 7 LLK = -1123.16 (3.763832)
## # 8 LLK = -1121.775 (1.385447)
## # 9 LLK = -1121.157 (0.617467)
## # 10 LLK = -1120.82 (0.3373408)
## # 11 LLK = -1120.592 (0.2281166)
## # 12 LLK = -1120.413 (0.178957)
## # 13 LLK = -1120.262 (0.1510623)
## # 14 LLK = -1120.13 (0.1311716)
## # 15 LLK = -1120.016 (0.11493)
## # 16 LLK = -1119.915 (0.1009953)
## # 17 LLK = -1119.826 (0.08892939)
df <- function(model) {
m <- length(model$mmu)
d <- length(model$mmu[[1]])
(m - 1) + m * d + m * d*(d+1) / 2
}
aic.gm <- function(model, YY) {
2 * (llk(model, YY) - df(model))
}
bic.gm <- function(model, YY) {
N <- nrow(YY)
2 * llk(model, YY) - df(model) * log(N)
}
Test df, AIC, BIC
model <- em.gauss(3, YY, do.plot=FALSE)
df(model)
## [1] 17
aic.gm(model, YY)
## [1] -2273.804
bic.gm(model, YY)
## [1] -2335.103
ic.mod.sel <- function(m.max, n.tries=10, do.plot=FALSE) {
aic <- numeric(m.max)
bic <- numeric(m.max)
for (m in 1:m.max) {
cat('m =', m, '\n')
model <- fit.em.gauss(m, YY, n.tries=n.tries, verbose=FALSE)
# model <- em.gauss(m, YY)
aic[m] <- aic.gm(model, YY)
bic[m] <- bic.gm(model, YY)
}
if (do.plot) {
matplot(cbind(aic, bic), type='b', pch=c(1,2), col=c('blue', 'red'), lty=c('solid'),
main='', xlab='m')
}
legend("topleft", bty = "n",
legend = c('AIC', 'BIC'),
col = c('blue', 'red'),
lty = c('solid', 'solid'))
return(list(aic=aic, bic=bic))
}
set.seed(1)
ic.mod.sel(7, 20, do.plot=TRUE)
## m = 1
## m = 2
## m = 3
## m = 4
## m = 5
## m = 6
## m = 7
## $aic
## [1] -2589.593 -2282.528 -2273.093 -2276.686 -2265.537 -2264.694 -2276.381
##
## $bic
## [1] -2607.623 -2322.192 -2334.392 -2359.619 -2370.105 -2390.897 -2424.218
plot.classification <- function(model, YY, plot.densities=TRUE) {
N <- nrow(YY)
m <- length(model$pp)
ZZ <- matrix(unlist(model$zz), N, m)
xlim <- c(min(YY[,1]), max(YY[,1]))
ylim <- c(min(YY[,2]), max(YY[,2]))
if (plot.densities) {
plot.mixture(model, YY, main=paste('m =', m))
} else {
plot(YY, type='n', main=paste('m =', m))
}
class.labels <- apply(ZZ, 1, which.max)
for (ell in 1:m) {
par(new=TRUE)
plot(YY[class.labels==ell,], col=(1+ell), xlab='', ylab='', xlim=xlim, ylim=ylim, axes=FALSE)
}
}
set.seed(1)
model <- fit.em.gauss(2, YY, 20, verbose=FALSE)
ZZ <- matrix(unlist(model$zz), 272, 2)
head(ZZ)
## [,1] [,2]
## [1,] 1.000000e+00 3.390575e-09
## [2,] 1.609446e-09 1.000000e+00
## [3,] 9.999897e-01 1.029217e-05
## [4,] 9.543642e-06 9.999905e-01
## [5,] 1.000000e+00 1.849258e-21
## [6,] 6.361447e-03 9.936386e-01
head(rowSums(ZZ))
## [1] 1 1 1 1 1 1
# matplot(ZZ)
plot.classification(model, YY)
par(mfrow=c(2,2))
for (ell in 2:5) {
set.seed(1)
plot.classification(fit.em.gauss(ell, YY, 30, verbose=FALSE), YY, plot.densities=FALSE)
}
par(mfcol=c(1,1))
Mclustlibrary('mclust')
## Package 'mclust' version 5.2.3
## Type 'citation("mclust")' for citing this R package in publications.
set.seed(1)
par(mfrow=c(2,2))
for (ell in 2:5) {
set.seed(1)
M <- Mclust(YY, G=ell, modelNames='VVV')
plot(M, what='classification')
}
par(mfcol=c(1,1))
par(mfrow=c(2,2))
for (ell in 2:5) {
set.seed(1)
M <- Mclust(YY, G=ell, modelNames='VVV')
plot(M, what='uncertainty')
}
par(mfcol=c(1,1))
par(mfrow=c(2,2))
for (ell in 2:5) {
set.seed(1)
M <- Mclust(YY, G=ell, modelNames='VVV')
plot(M, what='density', type = 'image', col = 'dodgerblue3', grid = 100)
# plot(M.2.VVV, what='density', type = 'persp')
}
par(mfcol=c(1,1))
set.seed(1)
model.2 <- fit.em.gauss(2, YY, n.tries=20, verbose=FALSE)
plot.mixture(model.2, YY)
df(model.2)
## [1] 11
# aic.gm(model, YY)
bic.gm(model.2, YY)
## [1] -2322.192
llk(model.2, YY)
## [1] -1130.264
M.2.VVV <- Mclust(YY, G=2, modelNames='VVV')
M.2.VVV$df
## [1] 11
M.2.VVV$bic
## [1] -2322.192
M.2.VVV$loglik
## [1] -1130.264
set.seed(1)
my.bic <- ic.mod.sel(7, 20, do.plot=TRUE)$bic
## m = 1
## m = 2
## m = 3
## m = 4
## m = 5
## m = 6
## m = 7
set.seed(1)
mc.bic <- mclustBIC(YY, modelNames='VVV')
cat('Our BICs:\n')
## Our BICs:
print(my.bic)
## [1] -2607.623 -2322.192 -2334.392 -2359.619 -2370.105 -2390.897 -2424.218
cat('Mclust\'s BICs:\n')
## Mclust's BICs:
print(as.numeric(mc.bic))
## [1] -2607.623 -2322.192 -2333.894 -2359.216 -2385.288 -2398.974 -2426.488
## [8] -2434.990 -2447.286
matplot(cbind(mc.bic[1:7], my.bic), pch=c(1,2), type='b')
MclustBIC <- mclustBIC(YY)
plot(BIC)
summary(BIC)
## Best BIC values:
## EEE,3 EEE,4 VVE,2
## BIC -2314.386 -2320.207047 -2320.432965
## BIC diff 0.000 -5.820699 -6.046616
BIC
## Bayesian Information Criterion (BIC):
## EII VII EEI VEI EVI VVI EEE
## 1 -4024.721 -4024.721 -3055.835 -3055.835 -3055.835 -3055.835 -2607.623
## 2 -3452.998 -3458.300 -2354.601 -2350.607 -2352.618 -2346.065 -2325.220
## 3 -3377.712 -3336.542 -2323.008 -2332.698 -2332.204 -2342.371 -2314.386
## 4 -3230.246 -3245.732 -2323.676 -2331.829 -2334.756 -2343.068 -2320.207
## 5 -3149.395 -3128.214 -2337.696 -2348.300 -2355.891 -2374.307 -2336.975
## 6 -3081.411 -3067.559 -2338.118 -2363.112 -2357.725 -2372.748 -2347.297
## 7 -2990.335 -2998.497 -2356.461 -2370.061 -2375.851 -2393.101 -2361.206
## 8 -2978.090 -2991.851 -2371.809 NA -2395.989 NA -2376.917
## 9 -2899.779 -2920.951 -2388.629 NA -2399.083 NA -2393.728
## EVE VEE VVE EEV VEV EVV VVV
## 1 -2607.623 -2607.623 -2607.623 -2607.623 -2607.623 -2607.623 -2607.623
## 2 -2324.273 -2322.972 -2320.433 -2329.116 -2325.416 -2327.598 -2322.192
## 3 -2322.690 -2322.094 -2332.433 -2338.986 -2329.352 -2335.618 -2333.894
## 4 -2340.044 -2332.885 -2354.890 -2336.750 -2342.472 -2344.720 -2359.216
## 5 -2354.462 -2352.661 -2376.045 -2356.225 -2366.193 -2374.132 -2385.288
## 6 -2368.534 -2364.464 -2392.133 -2371.732 -2387.445 -2387.218 -2398.974
## 7 -2381.415 -2375.217 -2398.950 -2392.963 -2384.183 -2413.721 -2426.488
## 8 -2401.464 NA NA -2385.815 -2404.950 -2434.531 -2434.990
## 9 -2404.102 NA NA -2418.305 -2428.351 -2431.113 -2447.286
##
## Top 3 models based on the BIC criterion:
## EEE,3 EEE,4 VVE,2
## -2314.386 -2320.207 -2320.433